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Abstract 

In this article we report the results of global structural optimization of the Si(114) sur- 
face, which is a stable high-index orientation of silicon. We use two independent procedures 
recently developed for the determination of surface reconstructions, the parallel-tempering 
Monte Carlo method and the genetic algorithm. These procedures, coupled with the use of 
a highly-optimized interatomic potential for silicon, lead to finding a set of possible models 
for Si(114), whose energies are recalculated with ab-initio density functional methods. The 
most stable structure obtained here without experimental input coincides with the struc- 
ture determined from scanning tunneling microscopy experiments and density functional 
calculations by Erwin, Baski and Whitman [Phys. Rev. Lett. 77, 687 (1996)]. 
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1 Introduction 



While the low-index semiconductor surfaces have dominated the interest of surface scientists 
for several decades, at present a considerable amount of work is being dedicated to high-index 
orientations. Since the high- index surfaces exhibit more diverse structural and electronic 
properties, the adept use of these properties can constitute the key for various technological 
applications, in particular for the fabrication of devices at length scales where lithographic 
techniques are not applicable. Controlling certain physical processes (e.g., growth of nanos- 
tructures in a preferred direction) depends in part on having knowledge of the structure of 
the substrate surface. The main technique for investigating atomic-scale features of surfaces 
is scanning tunelling microscopy (STM), although STM alone is only able to provide "a range 
of speculative structural models which are increasingly regarded as solved surface structures" 
PQ. A common procedure for finding the reconstructions of silicon surfaces consists in a 
combination of STM imaging and density functional calculations as follows. Starting from 
the bulk truncated surface and taking cues from the experimental data, one proposes several 
atomic models for the surface reconstructions. These models are then relaxed using electronic 
structure methods; at the end of relaxation, surface energies and STM images are computed 
for each structural model. A match with the experimental STM data is identified based on 
the relaxed lowest-energy structures and their simulated STM images (e.g., (21 El El)- As 
described, the procedure is heuristic, as one needs to rely heavily on physical intuition when 
proposing good candidates for the lowest energy reconstructions of high-index surfaces. 

Treating the reconstruction of semiconductor surfaces as a problem of global optimiza- 
tion, we have recently developed a parallel tempering Monte Carlo procedure for studying 
the structure and thermodynamics of crystal surfaces [5], as well as a genetic algorithm for 
structure determination The use of such methods can help avoid situations in which the 
actual physical reconstruction of a high-index surface is not part of the set of heuristic models 
that are considered for computation of surface energies and comparison with experimental 
data. Given that there are examples of semiconductor surfaces (e.g., [31 E]) for which the 
initially proposed models did not withstand further scientific scrutiny from different research 
groups, it appears worthwhile to perform searches for the structure of some of stable high- 
index surface orientations of silicon. One such surface is Si(114), reported to be as stable 
as the well studied low-index surfaces Si(OOl) and Si(lll) [Hj: given this stability of Si(114), 
it is somewhat surprising that this surface has not attracted more interest, at least from a 
technological perspective. There are few studies of Si(114) to date, which include the pioneer- 
ing study reporting on the atomic configuration 5 , and two recent works reporting on the 
electronic structure of this surface JHj • 

Based on scanning tunelling microscope (STM) images combined with density functional 
calculations, two atomic models [(2 x 1) and c(cx 2)], were proposed for the Si(114) orientation 
P3- These models have very similar bonding topology, differing only in terms of dimerization 
pattern of their surface. The surface energies of the two models are also similar, as both can 
be found on sufficiently large areas of the scanned samples To our knowledge, so far Ref. 
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[5] represents the only proposal for the structure of Si(114). The purpose of this article is to 
present several other model candidates for the structure of Si(114), models that are likely to 
be experimentally observed on this surface. Addressing the problem of atomic structure from 
a different perspective than the previous reports [31 E], we perform stochastic searches for the 
global minimum configuration of this surface. As we shall see, the lowest energy configuration 
(at zero Kelvin) obtained here from purely theoretical means is consistent with the original 
proposal [SJ; however, the global search methods provide several other structural models with 
low surface energy, which could be relevant in various experimental conditions. The remainder 
of this paper is organized as follows. Section II presents the two recently developed global 
search algorithms, the parallel tempering Monte Carlo methods and the genetic algorithm for 
structure determination. While brief descriptions of these methods are provided in Sec. II, 
we refer the reader to our recent works [51 E] for full details related to their implementation. 
The results of the structural optimization for Si(114) surface are presented and discussed in 
Section III, and our conclusions are outlined in the last section. 

2 Methods 

2.1 Parallel-tempering Monte Carlo method 

The reconstructions of semiconductor surfaces are determined not only by the efficient bonding 
of the surface atoms, but also by the stress created in the process (Hj. Therefore, we retain 
a large number of subsurface atoms when performing a global search for the lowest energy 
configuration: this way the surface stress is intrinsically considered when reconstructions are 
sorted out. The number of local minima of the potential energy is also large, as it scales 
roughly exponentially [HI El with the number of atoms involved in the reconstruction; by 
itself, such scaling requires the use of fast stochastic search methods. One such method is the 
parallel-tempering Monte Carlo (PTMC) algorithm |1UI lllj . which was shown to successfully 
find the reconstructions of a vicinal Si surface when coupled with an exponential cooling j^j. 
Before outlining the procedure, we discuss briefly the computational cell and the empirical 
potential used. 

The simulation cell (of dimensions 3a x ay/2 in the plane of the surface) has a single- face 
slab geometry with periodic boundary conditions applied in the plane of the surface, and no 
periodicity in the direction normal to it. The "hot" atoms from the top part of the slab (10- 
15 A thick) are allowed to move, while the bottom layers of atoms are kept fixed to simulate 
the underlying bulk crystal. The area of the simulation cell and the number of atoms in the cell 
are kept fixed during each simulation. Under these conditions, the problem of finding the most 
stable reconstruction reduces to the global minimization of the total potential energy V(x) of 
the atoms in the simulation cell (here x denotes the set of atomic positions). In terms of atomic 
interactions, we are constrained to use empirical potentials because the highly accurate ab- 
initio or tight-binding methods are prohibitive as far as the search itself is concerned. Since 
this work is aimed at finding the lowest energy reconstructions for arbitrary surfaces, the 
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choice of the empirical potential is important. After numerical experimentation with several 
empirical models, we chose to use the highly optimized empirical potential (HOEP) recently 
developed by Lenosky et al. ;12 . HOEP is fitted to a large database of ab-initio calculations 
using the force-matching method, and provides a good description of the energetics of all 
atomic coordinations up to Z = 12. 

The parallel tempering Monte Carlo method (also known as the replica-exchange Monte- 
Carlo method) consists in running parallel canonical simulations of many statistically inde- 
pendent replicas of the system, each at a different temperature T\ < T2 < . . . < T]y. The 
set of N temperatures {Tj, i = 1,2, N} is called a temperature schedule (or schedule for 
short). The probability distributions of the individual replicas are sampled with the Metropo- 
lis algorithm [E], although any other ergodic strategy can be employed ,14 3 . Irrespective of 
what sampling strategy is being used for each replica, the key feature of the parallel tempering 
method is that swaps between replicas of neighboring temperatures Tj and Tj (j = i± 1) are 
proposed and allowed with the conditional probability |ll)L ITT] given by 



min 




where V(xi) represents the energy of the replica i and is the Boltzmann constant. The 
conditional probability (^Q) ensures that the detailed balance condition is satisfied and that 
the equilibrium distributions are the Boltzmann ones for each temperature. 

In the limit of low temperatures, the PTMC procedure allows for a geometric temper- 
ature schedule ^3E3- ^o snow this, we note that when the temperature drops to zero, the 
system is well approximated by a multidimensional harmonic oscillator, so the acceptance 
probability for swaps attempted between two replicas with temperatures T < T' is given by 
the incomplete beta function law jTH] 

O rl/{l+R) 

= (2) 

where d denotes the number of degrees of freedom of the system, B is the Euler beta function, 
and R = T'/T. Since it depends only on the temperature ratio R, the acceptance probability 
P]) has the same value for any arbitrary replica running at a temperature Tj, provided that 
its neighboring upper temperature Tj + i is given by Tj + i = RT-i. The value of R is determined 
such that the acceptance probability given by Eq. © attains a prescribed value p. Thus, 
the (optimal) schedule that ensures a constant probability p for swaps between neighboring 
temperatures is a geometric progression: 

Ti = R l - x T min , l<i<N, (3) 

where T m j n = T\ is the minimum temperature of the schedule. 

The typical Monte Carlo simulation done in this work consists of two main parts that 
are equal in terms of computational effort. In the first stage of the computation, we perform 
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a parallel tempering run for a range of temperatures [T m j n , T max ]. The configurations of 
minimum energy are retained for each replica, and used as starting configurations for the 
second part of the simulation, in which replicas are cooled down exponentially until the 
largest temperature drops below a prescribed value. As a key feature of the procedure, the 
parallel tempering swaps are not turned off during the cooling steps. Thus, in the second 
part of the simulation we are in fact using a combination of parallel tempering and simulated 
annealing, rather than a simple cooling. At the k-th. cooling step, each temperature from the 
initial temperature schedule {Tj, i = 1,2,..., A?"} is decreased by a factor which is independent 
of the index i of the replica, T± = (X}.T^ . Because the parallel tempering swaps are not 
turned off, we require that at any cooling step k all N temperatures must be modified by the 
same factor at in order to preserve the original swap acceptance probabilities. We have used 
a cooling schedule of the form jS] 

if) = alf ~ 1} = a k - l Ti (k > 1), (4) 

where Tj = T^ 1 ' and a = 0.85 . 

The third and final part of the minimization procedure is a conjugate- gradient optimiza- 
tion of the last configurations attained by each replica. The relaxation is necessary because 
we aim to classify the reconstructions in a way that does not depend on temperature, so 
we compute the surface energy at zero Kelvin for the relaxed slabs i, i = 1,2,..., A. The 
surface energy 7 is defined as the excess energy (with respect to the ideal bulk configuration) 
introduced by the presence of the surface: 

7 = (E m - n m e b )/A (5) 

where E m is the potential energy of the n m atoms that are allowed to move, e b = — 4.6124eV 
is the bulk cohesion energy given by HOEP, and A is the surface area of the slab. 

At the end of the simulation, we analyze the energies of the relaxed replicas. Typical 
plots showing the surface energies of the structures retrieved by the PTMC replicas are shown 
in Fig. ^a), for different numbers of particles in the computational cell. To exhaust all the 
possibilities for the numbers of particles corresponding to the supercell dimensions of 3a x a\/2 ; 
we repeat the PTMC simulation for different values of n ranging from 208 to 220, and look 
for a periodic behavior of the lowest surface energy as a function of n. For the case of Si(114), 
this periodicity occurs at intervals of An = 4, as shown in Fig. ^b). Therefore, the (correct) 
number of atoms n at which the lowest surface energy is attained is n = 216, up to an 
integer multiple of An. As we shall show in the next section, the repetition of the simulation 
for different values of n in the simulation cell can be avoided within a genetic algorithm 
approach. 

2.2 Genetic Algorithm 

Like the previous method, the genetic algorithm also circumvents the intuitive process when 
proposing candidate models for a given high-index surface. An advantage of this algorithm 
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Figure 1: (a) Surface energies of the relaxed parallel tempering replicas i, (0 < i < 31) with 
total number of atoms n = 216 (solid circles), 215 (open triangles), 214 (open circles) and 213 
(solid triangles). For clarity, the range of plotted surface energies was limited from above at 
100 me V/A 2 . (b) Surface energy of the global minimum structure showing a periodic behavior 
as a function of n, with a period of An = 4; this finding helps narrowing down the set of 
values for n that need to be considered for determining the Si(114) reconstructions that have 
a 3a x a\/2 periodic cell. 
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over most of the previous methodologies used for structural optimization is that the number 
of atoms involved in the reconstruction, as well as their most favorable bonding topology, can 
be found within the same genetic search. The development of a genetic algorithm (GA) for 
surface structure determination was motivated by its successful application for the structural 
optimization of atomic clusters ^3 |20j . 

This search procedure is based on the idea of evolutionary approach in which the mem- 
bers of a generation (pool of models for the surface) mate with the goal of producing the 
best specimens, i.e. lowest energy reconstructions. "Generation zero" is a pool of p different 
structures obtained by randomizing the positions of the topmost atoms (thickness d) , and by 
subsequently relaxing the simulation slabs through a conjugate-gradient procedure. The evo- 
lution from a generation to the next one takes place by mating, which is achieved by subjecting 
two randomly picked structures from the pool to a certain operation (mating) 0:(A,B) — ►C. 
The mating operation O produces a child structure C from two parent configurations A and 
B, as follows. The topmost parts of the parent models A and B (thickness d) are separated 
from the underlying bulk and sectioned by an arbitrary plane perpendicular to the surface. 
The (upper part of the) child structure C is created by combining the part of A that lies 
to the left of the cutting plane and the part of slab B lying to the right of that plane: the 
assembly is placed on a thicker slab, and the resulting structure C is subsequently relaxed. 

A mechanism for the survival of the fittest is implemented as a defining feature of the 
genetic evolution. In each generation, a number of m mating operations are performed. The 
resulting m children are relaxed and considered for the possible inclusion in the pool based 
on their surface energy. If there exists at least one candidate in the pool that has a higher 
surface energy than that of the child considered, then the child structure is included in the 
pool. Upon inclusion of the child, the structure with the highest surface energy is discarded in 
order to preserve the total population p. As described, the algorithm favors the crowding of 
the ecology with identical metastable configurations, which slows down the evolution towards 
the global minimum. To avoid the duplication of members, we retain a new structure only if 
its surface energy differs by more than 5 when compared to the surface energy of any of the 
current members p of the pool. Relevant values for the parameters of the algorithm are given 
in 0: 10 < p < 40, m = 10, d = 5A , and 5 = 10- 5 meV/A 2 . 

We have developed two versions of the algorithm. In the first version, the number 
of atoms n is kept the same for every member of the pool by automatically rejecting child 
structures that have different numbers of atoms from their parents (mutants). In the second 
version of the algorithm, this restriction is not enforced, i.e. mutants are allowed to be part 
of the pool: in this case, the procedure is able to select the correct number of atoms for the 
ground state reconstruction without any increase over the computational effort required for 
one single constant-re run. The results of a variable-n run are shown in Fig. |2f a) which shows 
how the lowest energy and the average energy from a pool of p = 30 structures decreases 
as the genetic algorithm run proceeds. The plot in Fig. HJa) displays typical features of the 
evolutionary approach: the most unfavorable structures are eliminated from the pool rather 
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Figure 2: (a) Evolution of the lowest surface energy (solid line) and the average energy (dash 
line) for a pool of p = 30 structures during a genetic algorithm (GA) run with variable n 
(210 < n < 222). (b) Evolution of the number of atoms n that corresponds to the model 
with the lowest energy from the pool, during the same GA run. Note that the lowest energy 
structure of the pool spends most of its evolution in states with numbers of atoms that are 
compatible with the global minimum, i.e. n = 212 and n = 216. 
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fast (initial steep transient region of the graphs) and a longer time is taken for the algorithm 
to retrieve the most stable configuration. The lowest energy structure is retrieved in less than 
500 mating operations. The correct number of atoms [refer to Fig. HJb)] is retrieved much 
faster, within approximately 100 operations. It is worth noting that even during the transient 
period, the lowest-energy member of the pool spends most of its evolution in a state with a 
number of atoms (n = 212) that is compatible with the global minimum structure. 

The two independent algorithms (PTMC and GA) presented briefly in this section are 
able to retrieve a set of possible candidates for the lowest energy surface structure. We use 
both of the algorithms in this work in order to assess how robust their structure predictions 
are. As it turns out, the two methods not only find the same lowest energy structures for each 
value of the total number of atoms n, but also most of the other low-energy reconstructions - 
a finding that builds confidence in the quality of the configuration sampling performed here. 
Since the atomic interactions are modelled by an empirical potential ^2], it is desirable to 
check the relative stability of different model structures using higher-level calculations based 
on density functional theory; the details of these calculations are presented next. 

2.3 Density functional calculations 

Using the methodologies described above, we build a database of model structures that are 
sorted according to the surface energy given by the HOEP ^2] interaction model. Since 
the empirical potentials may not give a reliable energetic ordering when a large number of 
structures are considered, we recalculate the surface energies of the models in the database at 
the level of density functional theory. The calculations where performed with the plane-wave 
based PWscf package [U, using the Perdew-Zunger [22] exchange-correlation energy. The 
slab geometry and the computational parameters are similar to the ones reported in [Sj; given 
the increase in computational speed over the last eight years, we used thicker slabs and a 
different sampling of the Brillouin zone. The cutoff for the plane-wave energy was set to 12 
Ry, and the irreducible Brilloiun zone was sampled with 4 k points. The equilibrium bulk 
lattice constant was determined to be a =5.4lA, which was used for all the surface calculations 
in this work. The simulation cell has the single-face slab geometry, with 24 layers of Si, and 
a vacuum thickness of 12 A. The bottom three layers are kept fixed in order to simulate the 
underlying bulk geometry, and the lowest layer is passivated with hydrogen. The remaining 
Si layers are allowed to relax until the forces become smaller than 0.025 eV/A. 

The surface energy 7 for each reconstruction is determined indirectly, by first considering 
the surface energy 70 of an unrelaxed bulk truncated slab, then by calculating the difference 
A7 = 7 — 70 between the surface energy of the actual reconstruction and the surface energy 
of bulk truncated slab that has the bottom three layers fixed and hydrogenated. The energy 
of the bulk truncated surface, as computed from a two-faced slab with 24 layers, was found 
to be 70 = 143 meV/A 2 . This indirect procedure for calculating the surface energies at the 
DFT level was outlined, for instance, in Ref. |2j. 
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3 Structural models for Si(114) 



At the end of the global search procedures, we obtain a set of model structures which we 
sort by the number of atoms in the simulation cell and by their surface energy. Since the 
empirical potentials may not be fully transferable to different surface environments, we study 
not only the global minima given by the model for different values of n, but also most of the 
local minima that are within 15 meV/A 2 from the lowest energy configurations. After the 
global optimizations, the structures obtained are also relaxed by density functional theory 
(DFT) methods as described in Sec. 12.31 The results are summarized in Table Q which will 
be discussed next. 

3.1 Results 

Table 1 lists the density of dangling bonds (db per area), as well as the surface energies of 
several different models calculated using the HOEP potential and DFT. The configurations 
have been listed in increasing order of the surface energies computed with HOEP, as this is 
the actual outcome of the global optimum searches. For reasons of space, we limit the number 
of structures in the Table at most six for each value of the relevant numbers of atoms in 
the simulation cell. However, when performing DFT relaxations we consider more structures 
than the ones shown in the table because we expect changes in their energetic ordering at the 
DFT level. The inclusion of a larger number of structures helps avoid excessive reliance on 
the empirical potential |12j . which is mainly used a fast way to provide physically relevant 
reconstructions (i.e. where each atom at the surface has at most one dangling bond). 

Table H and Fig. ^b) suggest that the most unfavorable number of atoms in the sim- 
ulation cell is n = 215, both at the level of HOEP and at the level of DFT. Therefore, it is 
justifiable to focus our attention on the other three values of n (n = 216, 214 and 213), which 
yield considerably lower surface energies. For each of these numbers of atoms, we present four 
low energy structures (as given by DFT), which are shown in Figs. EH3 These structures are 
not necessarily the same as those enumerated in Table Q as they are chosen based on their 
DFT surface energies. Since the global optimization has not been performed at the DFT 
level, the reader could argue that the lowest energy structure obtained after the sorting of the 
DFT-relaxed models may not be the DFT global minimum. While we found that a thorough 
sampling for systems with ~ 200 atoms is impractical at the DFT level, we have performed 
DFT relaxations for most of the local minima given by HOEP. Therefore, given the rather 
large set of structural candidates with different topological features considered here, the pos- 
sibility of missing the actual reconstruction for Si (114) is much diminished in comparison with 
heuristic approaches. 

We will now describe in turn the surface models corresponding to n = 216, 214 and 
213. After the DFT relaxation, the lowest energy model that we found has turned out to 
be the same as the one proposed by Erwin and coworkers jS], perhaps with the exception 
of different relative tilting of the surface bonds. The model is shown in Fig. Eta), and it 
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n Bond counting HOEP DFT 

(db/3a 2 V2) (meV/A 2 ) (meV/A 2 ) 

216 8 81.66 89.48 

8 83.16 90.34 

8 83.31 91.29 

8 83.39 88.77 

8 83.64 94.68 

8 84.42 92.16 

215 8 91.61 97.53 

8 91.82 95.30 

8 92.00 94.20 

11 92.46 98.73 

214 6 86.95 95.17 

10 87.32 99.58 

10 87.39 98.47 

10 87.49 93.88 

10 88.26 95.18 

213 4 89.46 90.43 

6 89.76 94.01 

4 90.07 90.85 

6 91.73 94.66 

7 93.99 90.48 



Table 1: Surface energies of different reconstructions for the Si(114) surface, sorted by the 
number of atoms n in the 3a x ay/2 periodic cell. The second column shows the number 
of dangling bonds (counted for structures relaxed with HOEP) per unit area. The last two 
columns list the surface energies given by the HOEP interaction model |12j and by density 
functional calculations (DFT) [2^ with the parameters described in text. 
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is characterized by the presence of dimers, rebonded atoms and tetramers occurring in this 
order along the (positive) [221] direction. These features have been well studied [5JEI1> and 
we shall not insist on them here. The surface energy of the most stable model for Si(114)-(2) 
reconstruction is 7 = 88.77 meV/A 2 . Although this surface energy is somewhat different 
from the previously reported value of 85 meV/A 2 (Hj, the discrepancy between these absolute 
values can be attributed to the somewhat different computational parameters (slab thickness, 
number of k points) and/or different pseudopotentials. 

It is notable that a different succession of the above-mentioned atomic scale features 
is also characterized by a low surface energy: specifically, dimers, tetramers and rebonded 
atoms (in this order along [221]), as shown in Fig. EJc), give a surface energy which is only ~2 
meV/A 2 higher than that of the Erwin el al. model shown in|2ta). This surface energy gap 
is apparently large enough to allow for another configuration [see Fig. Etb)] with a surface 
energy that lies between the values corresponding to the first two models described above. 
As shown in Fig. Otb), this new model has two consecutive dimer rows followed by a row 
of rebonded atoms, and arrangement that gives rise to surface corrugations of 0.4-0.5 nm. 
Remarkably, this corrugated model [Fig. Otb)] is almost degenerate with the planar, (2 x 1) 
structure shown in Fig. |2ta). The last panel of Fig. |3] shows another planar model of Si(114), 
made of dimers, rebonded atoms and inverted tetramers, with the latter topological feature 
distinguishable as a seven-member ring when viewed along the [110] direction. 

Dimers, rebonded atoms and tetramers also occur on low-energy structural models with 
n = 214, as shown in Fig. ^ The most favorable structure with n = 214 that we found 
[depicted in Fig. Eta)] has a 5-coordinated subsurface atom and a 4-coordinated surface atom 
per unit cell. These topological features are determined by the bonding of a subsurface atom 
with one of the atoms of a tilted surface dimer; the corresponding surface energy is 90.09 
meV/A 2 . Other structures with n = 214 atoms [examples shown in Fig. Etb)-(d)] generally 
have higher energies than models with n = 216, (refer to Table most likely because the two 
missing atoms lead to pronounced strains in the surface bonds. 

The analysis of simulation slabs with n = 213 atoms reveals novel atomic scale features. 
Energetically favorable configurations with n = 213 [I^a) and (c)] show a 5-atom ring on the 
surface stabilized by a subsurface interstitial, a structural complex that was first encountered 
in the case of Si(113) surface |2j. Structures in Figs. IHfa) and (c) differ in terms of the 
succession of the topological features along the [221] direction, i.e. dimers, 5-member rings, 
rebonded atoms [St a)] as opposed to dimers, rebonded atoms and 5-member rings 0(c)] . The 
model in Fig. Eta) is degenerate with the one shown in|SJb), as their relative surface energy 
is much smaller than the 1-2 meV/A 2 expected accuracy of the relative surface energies 
determined here. The reconstruction Etb) is very similar to the lowest energy structure in 
Fig. Eta) (achievable with n = 212): the only different feature is the extra atom lying in 
between two rebonded atoms and sticking out of the surface [refer to Fig.Etb)]. Likewise, the 
model in Fig. Etd) can be obtained from structure Etb) by adding one atom per unit cell in 
such a way that it bridges the two atoms of a dimer on one side, and rebonds on the other 
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side. 



3.2 Discussion 

The data in Tabled shows clearly that the density of dangling bonds at the Si(114) surface 
is, in fact uncorrelated with the surface energy. The lowest number of dbs per area reported 
here is 4, and it corresponds to n = 213 and 7 = 90.43 meV/A 2 at the DFT level. The 
optimum structure [St a)] , however, has twice as many dangling bonds but the surface energy 
is smaller, 88.77 meV/A 2 . Furthermore, for the same number of atoms in the supercell 
(n = 216) and the same dangling bond density (8d6/3a 2 \/2), the different reconstructions 
obtained via global searches span an energy interval of at least 5 me V/A 2 . These findings 
constitute a clear example that the number of dangling bonds can not be used as a criterion 
for selecting model reconstructions for Si(114); we expect this conclusion to hold for many 
other high-index semiconductor surfaces as well. 

The HOEP surface energy and the DFT surface energy also show very little correlation, 
indicating that the transferability of the interaction model ^2j for Si (114) is not as good as, 
for instance, in the case of Si(001) and Si(105) [5]. The most that can be asked from this 
model potential [12, is that the observed reconstruction is amongst the lower lying energetic 
configurations -which, in this case it is. We have also tested the transferability of HOEP for 
the case of Si(113), and found that, although the ad-atom interstitial models are not the 
most stable structures, they are still retrieved by HOEP as local minima of the surface energy. 
We found that the low- index (but much more complex) Si(lll)-(7 x 7) reconstruction is also a 
local minimum of the HOEP interaction model, albeit with a very high surface energy. Other 
tests indicated that, while the transferability of HOEP to the Si(114) orientation is marginal 
in terms of sorting structural models by their surface energy, this potential performs much 
better than the more popular interaction models |17l I24j . which sometimes do not retrieve 
the correct reconstructions even as local minima. Therefore, HOEP is very useful as a way to 
find different local minimum configurations for further optimization at the level of electronic 
structure calculations. 

A practical issue that arises when carrying out the global searches for surface reconstruc- 
tions is the two-dimensional periodicity of the computational slab. In general, if a periodic 
surface pattern has been observed, then the lengths and directions of the surface unit vectors 
may be determined accurately through experimental means (e.g., STM or LEED analysis): in 
those cases, the periodic vectors of the simulation slab should simply be chosen the same as 
the ones found in experiments. When the surface does not have two-dimensional periodicity, 
or when experimental data is difficult to analyze, then one should systematically test compu- 
tational cells with periodic vectors that are integer multiples of the unit vectors of the bulk 
truncated surface, which are easily computed from knowledge of crystal structure and surface 
orientation. There is no preset criterion as to when the incremental testing of the size of the 
surface cell should be stopped -other than the limitation imposed by finite computational 
resources; nevertheless, this approach gives a systematic way of ranking the surface energies 
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Figure 3: Structural models (top and side views) of Si(114)-(2xl), with n = 216 atoms 
per unit cell after relaxation with density functional methods [21] . The surface energy 7 
computed from first-principles is indicated for each structure, along with the corresponding 
value (in parentheses) obtained using the empirical potential ^2]. The darker shade marks 
the undercoordinated atoms. 
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(a) 7 = 90.09 meV/A 2 (92.62) (b) y = 95.17 meV/A 2 (86.95) 




(c) 7 = 95.18 meV/A 2 (88.26) (d) y = 97.70 meV/A 2 (93.60) 



Figure 4: Structural models (top and side views) of Si(114)-(2xl), with n = 214 atoms 
per unit cell after relaxation with density functional methods |21| . The surface energy 7 
computed from first-principles is indicated for each structure, along with the corresponding 
value (in parentheses) obtained using the empirical potential |12j . The darker shade marks 
the undercoordinated atoms, while the overcoordinated atoms are shown in white. 
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(a) y = 90.43 meV/A 2 (89.46) (b) y = 90.48 meV/A 2 (93.99) 




(c) y = 92.56 meV/A 2 (92.28) (d) y = 94.66 meV/A 2 (91.73) 



Figure 5: Structural models (top and side views) of Si(114)-(2xl), with n = 213 atoms per 
unit cell after relaxation with density functional methods 21 . The surface energy 7 computed 
from first-principles is indicated for each structure, along with the corresponding value (in 
parentheses) obtained using the empirical potential |12j . The darker shade marks the under- 
coordinated atoms, while the white atoms are either four-coordinated or overcoordinated. 
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of slabs of different areas, and eventually finding the global minimum surface structure. 

Motivated by a previous finding that larger unit cells can lead to models with very 
low surface energies (see, for instance the example of Si(105) 013); we have also performed 
global minimum search using GA for slabs of dimensions 6a x a\/2, which correspond to a 
doubling of the unit cell in the [221] direction. The ground state structure at the HOEP 
level found in this case is still the corrugated model E^b) with a surface energy of 7 = 81.66 
meV/A 2 . As a low- lying configuration we again retrieve the original model (Hj with 7 = 83.39 
meV/A 2 . Furthermore, we also find a several models that have surface energy in between the 
two values, characterized by the presence of different 3a x ay/2 structures in the two halves 
of the 6a x ay/2 simulation cell. This finding suggests that [110]-oriented boundaries between 
different 3a x ay/2 models on Si(114) are not energetically very costly: this is consistent with 
the experimental reports of Erwin et al, who indeed found (2x1) and c(2 x 2) structures 
next to one another pj. 

4 Concluding remarks 

In conclusion, we have obtained and classified structural candidates for the Si(114) surface 
reconstructions using global optimization methods and density functional calculations. We 
have used both parallel-tempering Monte Carlo procedure coupled with an exponential cooling 
[£], as well as the genetic algorithm [Jj- Both of the methods are used in conjunction with the 
latest empirical potential for silicon which has a better transferability in comparison with 
more popular potentials ^3121]. We have built a large database of structures (reported, in 
part, in Table ^) which were further optimized at the DFT level. The lowest energy structure 
that we found (Fig. Eta)) after the DFT relaxation is the same as the one originally reported 
for Si(114)-(2 x 1) in 0. 

In addition, we have discovered several other types of structures (refer to Figs. EJb), 
Etc) and^a) andEfa)) that are separated (energetically) by 1-2 meV/A 2 from the lowest 
energy model 0- Given that the relative surface energies at the DFT level have an error of ±1 
meV/A 2 , and that the experiments of Erwin et al. [H] already identified two reconstructions 
([(2 x 1) and c(2 x 2)] whose surface energies are within 1-2 meV/A 2 from one another, it is 
conceivable that some of the models in Figs. EHS1 could also be found on the Si(114) surface. 
This prediction could be tested, e.g., by high-resolution transmission electron microscopy 
experiments such as the ones reported recently for the Si(5512) surface j2^j. Low-energy 
electron diffraction experiments, as well as more STM measurements could also shed light on 
whether there exist other structural models on a clean Si(114) surface than initially reported 
in Ref. 0. 
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